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Summary. We consider the nonlinear propagation of light in an optical fiber 
waveguide as modeled by the anharmonic Maxwell-Lorentz equations (AMLE). 
The waveguide is assumed to have an index of refraction which varies period- 
ically along its length. The wavelength of light is selected to be in resonance 
with the periodic structure (Bragg resonance). The AMLE system considered 
incorporates the effects of non-instantaneous response of the medium to the 
electromagnetic field (chromatic or material dispersion), the periodic structure 
(photonic band dispersion) and nonlinearity. We present a detailed discussion 
of the role of these effects individually and in concert. We derive the nonlinear 
coupled mode equations (NLCME) which govern the envelope of the coupled 
backward and forward components of the electromagnetic field. We prove the 
validity of the NLCME description and give explicit estimates for the deviation 
of the approximation given by NLCME from the exact dynamics, governed by 
AMLE. NLCME is known to have gap soliton states. A consequence of our re- 
sults is the existence of very long-lived gap soliton states of AMLE. We present 
numerical simulations which validate as well as illustrate the limits of the the- 
ory. Finally, we verify that the assumptions of our model apply to the parameter 
regimes explored in recent physical experiments in which gap solitons were ob- 
served. 
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1. Introduction 

There is a great deal of current interest in nonlinear optical phenomena in peri- 
odic structures. This interest has been fueled by advances in fabrication methods 
for periodic media and in their potential for use as components in all-optical com- 
munication systems and computers. The potential for applications is due to the 
rich variety of phenomena which result from the interactions of sufficiently in- 
tense (nonlinear) electro-magnetic fields with the underlying (linear) dispersion 
characteristics of the periodic structure [47] . The reason one may envision the 
use of nonlinear periodic structures in optical devices stems from the observa- 
tion that one can achieve very strong dispersion of a light pulse over very short 
distances by arranging the wavelength of light and period of the medium to be 
appropriately resonant. At sufficiently high intensities, one then expects a bal- 
ance between nonlinear and dispersive effects over short distances, thus giving 
rise to a rich class of phenomena in structures of small physical dimensions. 

This paper is motivated by experiments and theory on nonlinear wave propa- 
gation in one-dimensional periodic structures. Our goal is to validate the nonlin- 
ear coupled mode equations (NLCME), a model commonly used to describe this 
situation, and to clarify the roles played by the various physical mechanisms. The 
experiments involve the propagation of intense light in an optical fiber waveguide 
whose core has a periodically varying index of refraction along the length of the 
fiber, a fiber grating [25]. Experimentalists have observed the formation of gap 
solitons, solitary- wave-like localized structures whose time- frequency parameters 
lie in the photonic band-gap associated with the background periodic structure. 
These are of potential interest for use in all optical storage devices, since they 
can, in principle, travel at arbitrarily low speeds. Theoretical work on nonlinear 
propagation in periodic structures goes back to work of Winful ct. al. [44, 45], 
and Chen and Mills [8] . Explicit gap soliton solutions were derived in the context 
of a slowly varying envelope theory by Christodoulides and Joseph [9] and in a 
more general form by Aceves and Wabnitz [1]. For surveys on aspects related 
to this paper, see de Sterke and Sipe [11], Brown and Eggleton [7] and Kurizki 
et. al. [17]. Experiments demonstrating the existence of gap solitons have been 
performed by Eggleton, Slusher and collaborators [13, 14, 15], and by Broderick 
and his collaborators [6, 34]. In two and three dimensions, Akozbek and John 
have formally derived envelope equations and examined their solitary waves nu- 
merically [3]. 

In the remainder of this section we give brief overview of the underlying 
physics and modeling assumptions. We also introduce the analytical and numer- 
ical results developed in this article. 

Electromagnetic wave propagation in a dielectric medium is described by 
Maxwell's equations together with an appropriate constitutive relation describ- 
ing how electromagnetic waves interact with matter. An optical fiber has a high 
index core and a slightly lower index cladding. This index configuration confines 
rays to the core (total internal reflection) or, from the wave perspective, the 
index profile provides a potential well with a ground state (core mode) having 
most of its energy confined to the core. In the regimes which interest us here, to a 
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very good approximation the energy distribution has a fixed transverse structure 
given by the core mode and one may think of the transverse core mode ampli- 
tude as varying with time t and distance along the waveguide, z. In addition to 
this geometric constraint, we incorporate the following effects: 

(i) Non-instantaneous response of the medium to the field: The polarization, 
P, is related to the electric field, E, via an anharmonic Lorentz oscillator model. 

(ii) Periodicity of the medium: Spatial periodicity of the medium is built in by 
allowing the coefficients of the anharmonic Lorentz oscillator to vary periodically 
in space. A schematic of the physical system is shown in Figure 1. 

(iii) Nonlinear effects at appropriate intensities: The implied relation between 
P and E is such that regions of higher intensity \E\ 2 have higher refractive index. 
This is a so-called focusing (Kerr) nonlinearity. The localized region of higher 
intensity effectively creates an attractive potential well. 

The effects of non-instantaneous response and spatial periodicity each give 
rise to dispersion, the property that waves of different wavelengths travel at 
different speeds. The type of dispersion due to (i) is called chromatic or material 
dispersion and that arising due to effect (ii) is called photonic band dispersion. 
This results from interference effects arising from reflection and transmission in 
the periodic structure. 




Fig. 1. A schematic drawing of an optical fiber with periodic refractive index variation 



A model incorporating the above geometric constraints and physical effects 
is a variant of the Anharmonic Maxwell-Lorentz system [5, 36, 24], which incor- 
porates the spatial periodicity. This system is displayed in (6) and we shall refer 
to it below as AMLE. 1 

While in a bare (homogeneous, non-grated) optical fiber, light injected at one 
end of the waveguide will propagate with little back-scatter, significant back- 
scattering will occur in the presence of a periodic refractive index. This effect 

1 In the nonlinear optics literature, the relation between polarization P and electric field E is 
often taken to have the form: P = x(* — t)E(t)(It + . . . . The class of models we have 
chosen gives the same envelope equation, NLCME, in the scaling regime considered but has 
the added feature that it conserves energy. Energy estimates are central to our proof of the 
validity of NLCME: Theorem 1. 
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is most pronounced when the wave-length of light is roughly twice the grating 
period, 2d, the case of Bragg resonance. In this case there is strong coupling 
between backward and forward waves. We will assume that the variation of the 
index of refraction about its mean is small and is denoted by a parameter e. In 
terms of this parameter we consider the following scaling regime: 

• amplitude of the field <~ 0(^/e), and 

• initial spectral support of the pulse is concentrated in a wave number 
range of width 0(e) about ±/c.b = ±ir/d. 

Therefore, the spatial structure of the fields E and P may be viewed as functions 
of the form ^/eA(ex)e lkBX , where A(y) is a localized function of y. We shall refer 
to this as the slowly varying envelope approximation (SVEA). A schematic of 
this scaling ansatz is shown in Figure 2. 




Fig. 2. A schematic of a wave packet under the SVEA, with envelope with width e _1 , ampli- 
tude y/e, carrier wavelength 2d and a plot of the index of refraction variation with variation 
of size e and period d. 



Under these assumptions, nonlinear coupled mode equations (NLCME, see 
equations (10)) can be derived which govern the forward E + and backward E- 
propagating electric field wave envelopes on time scales of order e _1 . Thus, in 
this regime, the fine scale grating oscillations are effectively averaged and the 
original mathematical description in terms of a nonlinear partial differential 
equation with spatially periodic coefficients is replaced by a constant coefficient 
dispersive nonlinear partial differential equation. 
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The main results of this paper are as follows: 
Characterization of Phenomena: The formation of long-lived coherent struc- 
tures (gap solitons) is the result of a balance between the effects of dispersion 
and nonlinearity. 2 The energy in a wave packet, with frequency content local- 
ized about the Bragg resonant frequency, resolves into backward and forward 
propagating waves. If the field amplitude is appropriately large relative to the 
the amplitude of periodic variations in the medium, then wave energy does not 
disperse and is localized in space. Nonlinearity generates ever higher harmonics 
which is manifested in wave steepening and apparent carrier shock formation; see 
Sects. 3.4 and 3.6. However, in the presence of material dispersion due to finite 
time response, a stable balance between dispersion and nonlinearity is achieved; 
no shocks form and one has long-lived stable gap solitons (see Sect. 3.5, The- 
orem 1, and the corollary of Sect. 3.7). We also verify (Appendix A) that pa- 
rameter ranges corresponding to experiments are described by our model and 
theorem. 

Analytical — Theorem 1: We prove that solutions to the initial value prob- 
lem for AMLE with finite energy and nearly monochromatic initial conditions, 
as above, give rise to solutions which are well-approximated, on appropriate 
time scales (C(e -1 )) by a superposition of amplitude-modulated backward and 
forward propagating plane waves. The backward and forward amplitudes sat- 
isfy NLCME. Further, we estimate the deviation of the NLCME approxima- 
tion to the AMLE solution. An important class of solutions of NLCME are 
so-called gap solitons? These are spatially localized nonlinear bound states [1] 
which have been observed experimentally. Our results imply the existence of gap 
soliton wave packets for AMLE on time scales 0(e~ 1 ), see Sect. 6. Our method 
follows previous rigorous studies of the validity of solutions to envelope equa- 
tions in approximating oscillatory, nearly monochromatic solutions to evolution 
PDE's in one space dimension; see, for example, Kirrman-Schncider-Mielke [27] 
and Pierce- Wayne [37] . A related method was presented in the context of dis- 
sipative equations is given by van Harten [41]. Important extensions of such 
methods have been developed in Donnat-Rauch [12], Joly-Metivier-Rauch [23], 
Lannes [31], and Schochet [39, 40] where multidimensional and multiphase prob- 
lems are treated. These paperes do not treat the resonant interactions of multi- 
phase waves with periodic media, although the formulation of [12], e.g., can 
(easily) be extended to deal with the present one-space-dimensional problem. 
Here, however, we provide a self-contained, elementary treatment of the one- 
dimensional problem. We also wish to note a recent general paper of Babin and 
Figotin [4] on wave interactions in periodic media. 

Numerical Simulations: We numerically simulate AMLE and NLCME and 

2 In contrast to the case of bare fiber, for which this balance is achieved over lengths of optical 
fiber on the order of tens of kilometers, for the periodic structures and intensities used in the 
above cited experiments, this balance occurs on a length scale of centimeters; see Appendix 
B. 

3 Here, we adopt a common usage of the term soliton as referring to a nonlinear bound state 
solution or solitary wave. The term originally and still often refers more specifically to 
nonlinear bound states arising in completely integrable systems. 
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systematically compare their computer generated solutions with a view toward 
checking the error estimates of Theorem 1. Initial data appropriate for Theo- 
rem 1 consists of a slow modulation of a highly oscillatory wave. In the main 
numerical example presented in Sect. 7, we have simulated the AMLE evolution 
with gap soliton wave packet initial data. We take data with an envelope whose 
full width at half maximum (FWHM) measured in wavelengths of light, is about 
60, see Figure 7. In one run, the simulation of AMLE took about 2| hours, 
while the corresponding simulation of NLCME took only a few minutes on a 500 
mHz Pentium III computer running Linux. The advantage is even larger when 
wavepackets with more oscillations are investigated. For parameter regimes of 
physical interest, it is probably infcasible to simulate the full AMLE, while it is 
quite simple to simulate the NLCME. In the physical experiments [13], pulses 
on the the order of 30 ps FWHM are observed, O(10 4 ) wavelengths. 

In our numerical simulations we identify three time regimes. The first is the 
time scale on which the coherent structure evolves as a gap soliton plus fluctu- 
ations satisfying the estimates of Theorem 1. The second is a longer time scale 
on which the wave envelope predicted by NLCME gives a accurate prediction of 
where the field energy is, but due to phase drift, the norm estimates of the error 
in Theorem 1 fail to hold. The third is a regime on which the wave envelope 
begins to steepen asymmetrically and radiate energy, leading to a decay of the 
gap soliton. A description of this process would require the inclusion of higher 
order nonlinear wave-steepening and dispersive corrections to NLCME. 

The structure of this paper is as follows: 

In §2.1 we introduce the anharmonic Maxwell-Lorentz model (AMLE) and 
inn §2.2 we describe and display the nonlinear coupled mode equations, discuss 
their mathematical structure, and state our main theorem (Theorem 1) relat- 
ing solutions of AMLE to those of NLCME. In §3 we discuss the effect of the 
nonlincarity, periodic structure, and material dispersion and describe the phys- 
ical effects of including, excluding and variously combining these mathematical 
features of the system. In §4 we present a derivation of NLCME from AMLE 
using the method of multiple scales and in §5 we discuss existence and unique- 
ness results for AMLE and NLCME, some of which are needed in §6 where 
we prove Theorem 1. In §7 we report on numerical simulations and careful sys- 
tematic comparison of solutions to AMLE and NLCME. In §8 we present a 
short summary followed by a discussion of issues meriting further investigation. 
The appendices contain a discussion of nondimensionalization and physical pa- 
rameter magnitudes, and details of dimensionless values used in the numerical 
simulations. 



Notation and Conventions 

Throughout this paper we make use of following notation: 

The symbols C, Cj are used to represent generic constants whose dependence 
on parameters is specified when of concern. 
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For a vector- valued function f(z), the L p -norm is given by 

ii/1i p = (/E^^w 1/p • (!) 

Here and throughout, spatial (z) integrals are taken over — oo < z < oo. The 
space L p is then the space of all functions / such that ||/|| is finite. 
The L°° norm is given by 

1 1 /I I oo = m ax ( ess sup |/j («) | ) (2) 

with the space L°° thus defined as the set of all (essentially) bounded functions. 
The H s norms may be defined as 

ii/ii- =t\\Mi- < 3 > 
fe=i 

The space H s is the space of all functions / such that H/H^s is finite, that 
is, the space of functions such that the function and its first s derivatives are 
square-integrable. 

Finally, for a given Banach space X, with norm || • \ \x, we define 

C([0,T);X) 

to be the set of functions / : t \— > f(t) e X which are continuous for t € [0, T) 
with values in X. 



2. The AMLE and NLCME Equations 

2.1. AMLE, nondimensionalization and parameter regimes 

In this subsection we introduce AMLE with physical parameters and then in- 
troduce its nondimensional form. We then discuss parameter regimes associated 
with the above mentioned experiments. 

We take as our basic model a one-dimensional electromagnetic system sat- 
isfying Maxwell's equation, with the polarization governed by an anharmonic 
Lorentz oscillator model, 4 henceforth referred to as the anharmonic Maxwell- 

4 Our results and analysis apply to the generalization of this model where we take P to be a 
weighted sum of N polarizations, Pi, corresponding to different molecular excitation modes 
of the material: 

N 

P = Y J Pi ; u7 2 d 2 t Pi + (l - 2An t cos(2fe B2 )) P, - fcP? = eox^ E . 
i=l 

This model can be viewed as a nonlinear generalization of the Sellmeier model (f32); see [2], 
With this particular modeling of the polarization, AMLE has the important property of 
being an energy conserving system. This structure gives rise to energy estimates which 
are central to our proofs of well-poscdncss of AMLE and of the validity of NLCME as an 
approximating envelope equation; see Sects. 5 and 6. 
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Lorentz equations (AMLE) [5, 24, 36] 

u d t D = d z B; 8 t B = d z E ; (4a) 
D = e E + P ■ (4b) 

tio 2 d?P+ (l - 2Ancos(2k B z)j P- 4>P 3 = e X (1) E . (4c) 

Here, E is the electric field, B is the magnetic field, P is the polarization, and D 
is the electric displacement, eo and denote, respectively the permitivity and 
permeability of free space, and x^ is the linear polarizability of the medium. 
Recall that e ^o = c ~ 2 > where c is the vacuum speed of light. An measures the 
strength of the grating. We shall also write: 

An = ev , (5) 

where e measures the size of the index modulation and v is of order one and is 
introduced in order to make explicit how the spatially periodic structure rears 
its head in the envelope approximation, NLCME, to be derived below. 
The spatial period of the medium is d and is expressed in terms of 

d 

Since we are interested in the propagation of light with wavelength equal to the 
Bragg wavelength, we set 

A = 2d , 

where d denotes the period of the grating. 

In Appendix A, we eliminate the magnetic field B from this system, and 
nondimensionalize these equations. There, nondimensional dependent variables 
are primed, but here we drop primes for simplicity of notation. From (127), 
(125b), and (125c), we obtain: 

8 2 t D - d 2 z E ; (6a) 
D = E + P ; (6b) 
uio 2 d?P + (1 - 2ev cos(2k B z))P - ^P 3 = (n 2 - l)E . (6c) 

uiq is a dimensionless frequency and 4>, a dimensionless measure of the degree of 
the nonlinearity. The limit of instantaneous polarization is achieved by taking 
the parameter uj — > oo. This gives the relation 

P = P(E) = (1 + 2ev cos(2k B z))E + X {3) E 3 + ... 

with nonlinear polarizability x^ an d, in this case, equation (6) reduces to the 
scalar nonlinear wave equation. 



d 2 {n 2 E + 2evcos(2k B z)E + X i3) E 3 ) = d 2 E 



(7) 
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2.2. NLCME and Main Results 

The nonlinear coupled mode equations are introduced by considering slow mod- 
ulations to solutions of the anharmonic oscillator model in which the photonic 
structure and nonlinearity are ignored. When e = and 4> = system (6) sup- 
ports plane wave solutions of the form E = E±e l ^ ±kz ~ u ' t \ where k = k{uj) is 
the dispersion relation (see Sect. 3) and E± are constants. A similar statement 
holds for P. In the scaling regime described in the introduction, in which non- 
linear effects and spatial periodicity are allowed, and where the carrier wave 
has wavenumber Ub and frequency lub = w(fcs), we seek coupled and slowly 
modulated backward and forward plane wave solutions of the form 

e nlcme\ ^ y/i (E+(Z,T)e i( - kBZ - WBt) + E_(Z,T)e- i( - kBZ+WB ^ + c.c) ( 1 

P NLCME J V ' \7B. 

(8) 



where c.c. denotes the complex conjugate of the previous expression and 75 = 
7(w_b) is a constant. Here, Z and T are "slow variables:" 

Z = ez, T = et , (9) 

and E + and E- satisfy equations of the form 5 

i(d T E + +v g dzE + ) + K E-+r(\E + \ 2 + 2\E_\ 2 )E+=0 ; (10a) 
t(d T E^-v g d z E-) + KE + + r(\E-\ 2 + 2\E+\ 2 )E_ =0 . (10b) 

Here, k is a coupling parameter (proportional to v induced by the grating), v g 
is the group velocity of the linear dispersive wave at frequency wg, and r is the 
nonlinear coupling parameter (proportional to cf>). The explicit expressions for 
these coefficients are displayed in Sect. 4, during the derivation of (10); see (61). 

The expression in (8) for E £ and P £ is a formal approximate solution to 
AMLE satisfying the "nearly monochromatic" initial condition: 

pj^ = 0)) =Vi(Eo + (ez)€ + e ik ' + E -(ez)€-e- ik '+cc)+0( E ) (11) 

where v± are constant two-component vectors. 
We prove the following result in Sect. 6: 



Beginning with a three dimensional Maxwcll-Lorcntz model in fiber geometry, it is possible 
to derive similar nonlinear coupled mode equations, with one difference being that the term 
r(|B ± | 2 + 2\E^\ 2 )E± is replaced by one of the form (r s |_B ± | 2 + 2F X \E^\ 2 )E±, where T s 
and r x are the nonlinear self-phase modulation and cross-phase modulation coefficients, 
and depend on certain integrals of the transverse modes of the waveguide. [1, 35] Another 
difference one finds is that the transverse potential, defined by the refractive index profile, 
induces waveguide or modal dispersion. Thus, a more complete description of the physics 
leads to corrections to the free space dispersion relation due to material dispersion, photonic 
band dispersion and modal dispersion, as well. The multidimensional analyses of [12, 23, 31, 
40, 39], which assume "almost" plane wave solutions, do not appear to generalize (easily) to 
the waveguide problem. 
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Theorem 1 Consider AMLE with a general nonlinearity satisfying Hypothe- 
sis 2 of Sect. 6. There exists eo > such that for any T > and any < e < £o, 
the solution ( E pi MLE ) of AMLE with data (11) belonging to H 3 is well approxi- 

A M LE 

mated by a solution of NLCME in the sense that for allt E [0,T /e] the following 
estimate holds: 



'E 



AMLE 



\ AMLE / 



Enlcme\ 

^ NLCME J 



< C(T ; w , v, n)e 



(12) 



H 1 



We note that due to Sobolev's inequality, |/(a:)| < C||/||#i, a small error in 
the H 1 norm ensures a small pointwise error, so that the above statement gives 
uniform bounds on the error. 



3. AMLE, NLCME and Physical Phenomena 



The physical phenomena modeled by AMLE and NLCME result from compe- 
tition among: (i) nonlinear effects, (ii) dispersion due to finite time response of 
the medium to the field and (iii) dispersion due to reflection and transmission 
in a spatially periodic medium. This section is divided into subsections in which 
we study, by considering various choices of loq, e and <ft in (6), the action of these 
effects (terms in the equations) individually and in concert. 



3.1. Linear spatially homogeneous structure with instantaneous response: 

In this case, £ = = and uio = oo. Therefore P = (n 2 — l)E and the evolution 
is described by the classical wave equation: 

n 2 d?E = d 2 z E , (13) 

whose solutions are of the form 

E(z,t) = e + (z-t/n) + e-(z + t/n) , (14) 

corresponding to a a superposition of left and right moving waves which propa- 
gate without distortion. 

Alternatively, we can first seek elementary plane wave solutions E(z, t; k) — 
e -iuit+ikx _ ^ e ^en find that uu and k are related by the simple dispersion relation 
oj(k) = ±K Since the phase velocity, uj(k) /fc, is independent of k, all wavelengths 
travel at the same speed and we refer to (13) as nondispersive. Standard Fourier 
superposition of these plane waves yields the general solution (14). 
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3.2. Linear and homogeneous medium with finite time response 
In this case we have 

d 2 (E + P)=d 2 E; (15a) 
Lu^ 2 d 2 P + P=(n 2 - 1)E . (15b) 

We may still find plane wave solutions 

(£) - Q e *""~"" • (16) 

where k and u are related by the dispersion relation 

2 

n" 



1-1^- 



and 

n 2 -l 



k = <-o 7 / 2 , (17) 



1 - I — 

Mo 



2 • 



(18) 



In the relevant parameter regimes loq > u 2 and n 2 — 1 > 0. So, 7 is positive, 
corresponding to polarization in phase with the electric field. Note that in the 
Wo — > 00 limit, we recover the wave equation dispersion relation k = ±luu. 

A general solution may be constructed by superposition of these plane waves 
using the Fourier Transform. Using the method of stationary phase [43] , one can 
show that for initial data whose Fourier transform decays sufficiently rapidly, 
the amplitude of the solution decays as t~2 . 



3.3. Linear periodic structure, instantaneous response 

In this case, we take <f> = and £ / 0. We consider the case of instantaneous 
response, uj — 00, though the methods apply to finite wo as well. In this case 
we have the scalar one-dimensional wave equation with spatially periodic wave 
speed: 

(n 2 + 2e(n 2 - l)i/cos(2fc B z)) d 2 E = d 2 E . (19) 

In analogy with the scalar and spatially homogeneous wave equation, we seek 
solutions of the form: E(z, t) — e~ luJt Lp(z). This yields the Mathieu equation: 



-d 2 f{z) = uj 2 (n 2 + 2e(n 2 - l)j/cos(2fc B z)) ip(z) 



(20) 
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We now seek solutions of (20) of the form 



tp(z) = e iKz ^{z;K) 



K e [0, 2tt) 



(21) 
(22) 



where ip has the same periodicity as the medium. Therefore, ip(z; K) satisfies 
the boundary value problem: 



- (d z + iK) 2 i/>(z; K) = u? (n 2 + 2e{n 2 - l)u cos(2k B z)) ip(z; K) ; (23) 



For each K, there is a discrete set of eigen-solutions {ip m (z; K) : m — 1, 2, . . . } 
with corresponding eigenvalues {oj m (K) 2 : m = 1, 2, . . . }. As K varies over the 
interval [0, 27r), the functions ui^K) sweep out spectral (photonic) bands. These 
bands are separated by spectral (photonic band) gaps. The solutions 



E m (z, t; K) = e -^W+«^ m(z . K)j K e [0> 2 ^ m = i ; 2 , . . . . (25) 



are generalizations of plane wave solutions of the previous subsections. In con- 
trast to the homogeneous medium case, where the allowable set of frequencies 
varies over the entire real line, in the periodic case the allowable set of frequencies 
varies over selected bands. The band dispersion relations, {ui m {K) : m = 1, 2, ..}, 
play the role of the dispersion relation, oj(K), for the homogeneous medium 
(constant coefficient partial differential equation). Since the phase velocities 
uj, m {K) / K are not independent of K, we see that wave propagation in periodic 
media is dispersive. A generalization of Fourier superposition holds, enabling one 
to construct the general solution to the initial value problem for (19). A careful 
stationary phase analysis of this superposition formula can be made, yielding 
results on the spreading and temporal decay of solutions [30]. 

Thus, Floquet-Bloch theory gives a complete characterization of wave prop- 
agation in a linear periodic medium. However, the key to understanding the 
detailed properties of this propagation is a detailed knowledge of the band dis- 
persion functions, uj m (K). This is a difficult problem, in general. In the case 
when the periodicity is given by a small oscillation about its mean (e small), 
coupled mode theory [29] can be used to approximate the Floquet-Bloch spec- 
tral theory. This provides a satisfactory description of the wave propagation for 
large, but finite, times. We illustrate this for equation (19). The idea is that for 
e small, the solutions E m (z,t;K) should be well-approximated by plane waves 
of the unperturbed e — problem. Thus, we seek solutions of (19) in the form: 



and derive equations for the slowly varying functions E + (Z,T), E_(Z,T), en- 
suring that (26) is a good approximation of a solution for times, t, of order e -1 . 
This approximation and error estimates are derived systematically in Sects. 4 



1>(z + d;K) = $(z;K), d= — . 



(24) 
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and 6 in the nonlinear context of AMLE. In this linear setting the equations 
reduce to the linear coupled mode equations (cf. equations (10)): 

i (d T E+ + v g d z E+) + kE- = ; (27a) 
i {d T E_ - v g dzE-) + kE+ = , (27b) 



where v g = u>' = n 1 is the group velocity (which happens to agree here with 

k{n 2 -l)u 



the phase velocity, u)(k)/k), and k 



2n 

The opening of the first "photonic band gap" can be deduced from (27). 
Seeking solutions to (27) of the form 

E +^j = e i(QZ-a(Q)T) ^ ^ 2g ^ 

we find 

Q 2 (Q) = n- 2 Q 2 + n 2 . (29) 

The photonic band gap is pictured in Figure 3 which clearly shows a region 
of excluded frequencies centered around J? = 0. For J? in the gap, Q is imag- 
inary, indicating that frequencies in Bragg resonance with the grating cannot 
propagate. 




Fig. 3. The dispersion relation for the linearized coupled mode equations, showing the spectral 
gap. 



Finally, combining (28) and (29) with (26) gives the following approximation 
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to a band edge Floquet-Bloch generalized plane wave: 
E(z,t; K) \ K =k B +sQ 

= £ +e i[(bB-eQ)z-(u> B +en(Q))t] + g_ e -i[(k B +eQ)z+(u B +en(Q))t] + _ (39) 



3.4- Nonlinearity, instantaneous response and no periodic structure 
Here, we take uu — ► oo and e = in (6). The equations then reduce to: 

d 2 t (E + P)=d 2 z E ; (31a) 

and, for small E, 

P = P{E) = (n 2 -l)E + x {3) E 3 + ... . (31b) 

These may be combined to give 

d 2 D(E) = d 2 E (32a) 
D{E) = n 2 E + x (3) E 3 + ... . (32b) 

To study this system we first rewrite it in a more standard form. Let v = 
/C" 1 ^) = D(E). Then, (32) becomes d 2 v = d 2 z K{v). Introducing d z u = v, we 
obtain after one integration 

d 2 u = d z K{d z u) . (33) 

Equation (33) has the form of the equation governing the vibrations of a non- 
linear string, where the electric displacement, D(E), plays the role of the strain, 
d z u. 6 

A classical result of Lax [32] states that systems which are genuinely non- 
linear in the sense that JC"(0) ^ 0, or cquivalently D"(0) ^ 0, will develop 
singularities in finite time. Since D"(0) = 0, the quasilinear (32) does not satisfy 
the genuine nonlinearity condition, although D"'(0) ^ 0. Klainerman and Ma- 
jda [28] generalized Lax's result; if dg +1 ^D(0) ^ and the initial data is of size 
e then singularity formation takes place within a time interval of length (D(e~ p ). 

In particular, it follows from this result that for initial data of size 0{y/e) 
(see (8)), u(z,t) develops a singularity in its second derivatives within 0(e^ 1 ) 
time. Thus, E remains bounded but d z E tends to infinity at the singularity 
time. This is a shock type singularity. Specifically the results of [28] imply the 
following: 

6 The classical relation between tension, r and strain d z u, is derived via 



K(d z u) = d z u(l + (d z u) 2 y^r{d z u). 
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Theorem 2 Consider the quasilinear wave equation (32) with smooth initial 
data E(z,t = 0), d t E(z,t = 0) which are of order y/e. Then, there exists a 
finite and positive time, T(e) < Ce^ 1 such that 

sup ||£(-,t)||oo < oo , (34) 

0<t<T(e) 

while 

lim WdzEi-^Wcc + \\d t E(;t)\\°o = oo . (35) 

%/ 1 {e) 

Such carrier shock formation in the context of nonlinear optics has been discussed 
by heuristic arguments in [16, 19]. 

Figure 4 shows a simulation of the shocking process on a single carrier wave. 7 
Computational details are given in Sect. 7, and computational parameters for 
this and all subsequent numerical results are given in Appendix C. 




1 2 3 4 5 6 



X 

Fig. 4. Evolution from sinusoidal initial conditions (solid line) to near shock formation (dashed 
line, with shock location at the x) in Maxwell's equation with instantaneous nonlinear polar- 
ization. 



7 This computation is actually performed on the AMLE system with a very fast material 
response uo>l; sec Appendix C. The simulation is run to the time the shock "would have 
formed" in the absence of material dispersion. 
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3.5. Nonlinearity, finite time response and no periodic structure: 

The mathematical model in this case is AMLE, (6) with e = 0. Joly, Metivier, 
and Rauch [24] proved that the initial value problem for the full three dimen- 
sional AMLE, for some class of nonlinearities, does not develop singularities 
in finite time. In Sect. 5 we outline a proof of this result for our simpler one- 
dimensional model. Therefore, material dispersion, resulting from the finite re- 
sponse time (u>o < oo), inhibits shock formation by providing a mechanism for 
expelling high frequency modes away from the steepening regions. 

Numerical experiments suggest an interesting small dispersion limit as loq 
tends to infinity in (6c). Note that for < uo < oo the system is semilinear, but 
the limiting system is quasi-linear. 

Two computations with increasing values of u>o are shown in Figures 5 and 6 
at a short time after the shock formation time in the dispersionless (wo = oo) 
limit. As luq — > oo the number of oscillations increases and in some weak sense, 
the solution more closely approximates a weak solution to the Maxwell system 
with instantaneous nonlinear polarization. 




Fig. 5. Solution of AMLE after the "shock time" with small ujq. 



The small material dispersion (uj 3> 1) limit of AMLE is analogous to the 
small dispersion limit of the Korteweg-de Vries equation (KdV) , 

d t u + ud z u + ed^u = , (36) 
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0.25 




Fig. 6. Solution of AMLE after the "shock time" with ujq twice that of Figure 5. 



the equation of the free surface of an air/water interface in the regime of long 
waves of small amplitude. The dispersionlcss e — equation is the inviscid Burg- 
ers equation and is easily seen to develop shocks (singularities in derivatives of u) 
in finite time [43] . For e ^ solutions of KdV do not develop singularities [26] . 
For initial data which give rise to shock formation for e = 0, one observes, for e 
small, a scenario analogous to what we observe for AMLE in the limit of uj large. 
KdV is an integrable Hamiltonian system which is exactly solvable using the in- 
verse scattering transform (1ST) [18]. 1ST was used by Lax and Levcrmorc [33] 
and by Venakides [42] to study this small dispersion limit. In particular, the 
generation of oscillations is related to the dynamics of solitons. As in the case 
of KdV, for AMLE one observes the generation of solitary wave like oscillations 
as a result of carrier wave steepening. Computer simulations indicate that these 
solitary waves interact more strongly and generate radiation, a manifestation of 
AMLE's apparent nonintcgrability. 

3.6. Periodic structure with nonlinearity, instantaneous response 

In this case, the electric field is governed by 

8?D(E,z) = %E ; (37a) 
D(E, z) = (n 2 + 2svcos(2k B z)) E + c/}E 3 . (37b) 
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The multiple scales approach implemented in Sect. 4 formally yields an ex- 
pansion of the form (compare with (8)) 

E = V^J2 (^ m) (Z,T)e mlfes ( z -*/ n ) + E {m \Z : T)e- lmk ^ z+t l n ^ + c.c. + e*Ex , 



where — E±*\z,T) = E < ^ l \ez,et), m > 1 satishes a coupled system of 

infinitely many partial differential equations. This is in contrast to the case wo < 
oo, where the expansion is replaced by (8) involving only the two amplitudes E±^ 
at leading (0(y/ej) order. 

The reason for this difference can be seen by examining the equation for the 
correction, E\, which takes the form: 

(n 2 d 2 t -dl)E 1 = Y J [ A+(T, £) e <9(**-<"(fc)t) + A~{T, Z)eM kz+u}< - k ^ ] 

q ~ l (39) 

= YyA+{T,zy qk{z - t ' n ) + A-{T,zy qk{z+t ' n ^ ] . 

q>l 

The coefficients, A^ involve the unknown amplitudes E^ and their derivatives. 
In order for eE\ to be smaller than the first term in the expansion of E, it is 
necessary to remove all resonances from the right hand side. Resonances are 
excited by components of the right hand side which are plane waves of the 
homogeneous problem. If < uj n < oo, the unperturbed dynamics are dispersive 
(ui(qk) 7^ quj(k)). Therefore, the contributions to the above sum for q > 2 are 
nonresonant and the nonresonance condition implies coupled equations for E^ 

and E^ alone. In the case of instantaneous response (loq = oo: absence of 
material dispersion) all terms in the sum are resonant. Therefore, in order to 
preclude secular growth, we require = 0, m > 1. This yields a coupled 
system of infinitely many equations governing the evolution of the backward 
and forward amplitudes: E^\z,T), m > 1. We do not address the question 
of whether the approximate solution generated, via (38), is a convergent series 
which represents an approximate solution of Maxwell's equation. 

Indeed, the contrast we find between the dispersive (o>o positive and finite) 
and nondispcrsivc (ujq = oo) cases is consistent with the observations of the 
previous section concerning shock formation and therefore the generation of 
high frequency harmonics. 8 

Consequently, the NLCME system does not describe the evolution of the 
wave packet envelope for the system without material dispersion. Although the 
ratio of the effects of photonic band dispersion to material dispersion in the 

8 Note, however, that the simulations described in the previous section are not for wave- 
packet initial conditions. It is reasonable to ask whether the dispersion which comes from 
the photonic band structure is sufficient to regularize shocks. Preliminary direct simulations 
for the nondispersive (ojq = oo) limit employing numerical schemes designed to capture 
shock-like structures indicate that shocks very likely form in the carrier, though the envelope 
appears to evolve smoothly (E. Kirr, in progress). 



(38) 
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experiments of Egglcton, Slushcr et. al. [13, 14, 15] is of order 10 6 , we argue that 
nonlinearity rapidly (on a time scale of order loq 1 ) generates frequency content 
for which material dispersion is significant; see the appendix. As noted in 3.5, 
material dispersion regularizes the wave steepening by propagating nonlinearly 
generated frequencies, which are nearly resonant, away from a steepening front. 



3. 7. Periodic structure, nonlinearity and finite time response 

In this case, we have the full AMLE equations (6). We show in Theorem 1 that 
for small amplitude waves in the SVEA regime, solutions to AMLE are well 
approximated by solutions to the Nonlinear Coupled Mode Equations (NLCME). 

The NLCME have a well-known class of solutions known as "gap solitons" 
which are able to propagate through the fibers at any velocity between zero and 
the speed of light. We present them in the general form as derived by Aceves 
and Wabnitz [1]. The solutions depend on two parameters, \v\ < 1 and 5: 



E+ = sue 111 
E_ = -ae vn 



2r 



(sin 8) e lsa sech (6 - isS/2) 



2r 



A(smS) e lsa sech (0 + isS/2) 



(40a) 
(40b) 



where: 



VT^2 - 

6 = jK(smS){v- 1 Z -vT) 
s = sign (kT) ; 

e vn = ( — 



A 



a = 



1 



1 + v / 
jk{cosS)(v- 1 vZ -T) 



2(1 -v 2 ) 



,,20 



+ e 



-isS 



^29 _j_ gis& 



Combining this family of exact solutions to NLCME with Theorem 1, we 
have the following corollary: 

Corollary 1 The gap solitons approximate to 0(e) a family of long-lived solu- 
tions to the AMLE system for times of 0(e^ 1 ). 

The gap solitons solutions bear a striking resemblance to solitons of the 
Nonlinear Schrodinger equation (NLS). In fact, in the limit 8 -C 1, v <C 1, we 
may show that the gap soliton may be written as a normal mode of linear coupled 
mode equations, slowly modulated by an NLS soliton. To see this, we study the 
NLCME themselves under the SVEA limit. We assume 8 small, and look for 
approximate solutions to (10) of the form: 



E_ 



8A(8Z, ST, s 2 T)Ve^ QZ - QT ^ 



(41) 
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where V e i(QZ-nT) solveg the i inearize d NLCME. Then Q± = ±^Jk 2 + v 2 Q 2 
and A solves 



with 



ldT A+^!ld 2 c A + N\A\ 2 A (42) 



( = 5(Z-f2'(Q)T) , t = 5 2 T , 

r ( v 2 a Q 2 \ 



2 V" k 2 +v 2 Q 2 



NLS has spatially localized standing wave solutions of the form 



(43) 

and if we let Q = and A = | in this formula, then we recover exactly the 
leading term in the expansion of the gap soliton for v = and 5 -C 1. De Sterke 
and Sipe show additionally that the first two terms in the expansion of the gap 
soliton for small <5 and v correspond to the 0{5) and 0(S 2 ) terms in the multiple 
scales construction of solutions to AMLE [11] with small wavenumber Q. 

Therefore, we expect the following relationship among AMLE, NLCME and 
NLS. For 8 < S sufficiently small NLCME has a solution of the type (41), 
where ^4(C, t) satisfies NLS. The validity of NLS as an approximation to NLCME 
could be shown using the methods presented in Sect. 6 and [27]. This solution 
of NLCME, generated by NLS, gives rise to a solution of AMLE of the type (8), 
provided e < s(So) is sufficiently small. 

4. Derivation of the Nonlinear Coupled Mode Equations 

In this section we use the method of multiple scales [43] to derive the nonlinear 
coupled mode equations. We begin with the equation: 

d 2 (E + P)=d 2 E . (44a) 

We also specify a more general form for the nonlinear response in modeling the 
polarization in (6c): 

uj- 2 d?P+ (1 - 2evcos(2k B z))P + g(P,z) = (n 2 - l)E , (44b) 

where, for small values of P, 

g(P, z) — —<j)P 3 + higher order terms. (45) 

We expand the dependent variables in powers of e: 

E = e^E + e^E 1 +e r iE 2 + ... (46) 

P = e 1 iP {) + eip l +sip 2 + ... (47) 
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and expand the derivatives in terms of slow scales Z = ez and T = et: 

dt^d t + ed T and d x -» d z + ed z . (48) 

To derive the NLCME, it will be sufficient for us to consider the equations for 
the first two terms in the expansion, which may be written 



PoJ \2 V cos{2kBz)Po + (t>P$, , ' 



where 



and 



(49) 
(50) 



We now seek solutions order by order. 



0(e 1 ^ 2 ) At this order, the solution to the linear problem is given by (16), where 
fc and ui satisfy dispersion relation (17). As fc will be determined by the length 
scale of the Bragg grating structure, we prefer to solve (17) for w as a function 
of k: 

u 2 = l - (n 2 cj 2 + fc 2 ) ± ^{n 2 u 2 + fc 2 ) 2 - 4w 2 fc 2 . (53) 

Equation (53) has two roots corresponding to each choice of sign. In the limit as 
ojq — > oo, the root corresponding to the plus sign diverges to oo, while the root 
corresponding to the minus sign approaches the finite value n 2 u 2 — fc 2 , as noted 
in Sect. 3.1. We examine a pair of backward and forward propagating modes in 
Bragg resonance with the fiber and having slowly varying amplitudes: 

(p°) = (E+(Z,T)e i V CBZ -" Bt) +E-(Z,T)e- l(kBZ+ ^ Bt) + c.c.) Q ^ , (54) 

where 

ks = \ , (55) 
n 2 - 1 

1B= x _ ( ^ )2 , (56) 



and wb is a root of (53) for the minus sign choice. 



22 



R.H. Goodman, M.I. Weinstein, P.J. Holmes 



0(e 3 / 2 ) The equation at this order is 

£ °0 = + (2,cos(2fc B !)P o + 0Po 3 1 • 



Substituting in the solution to the 0(e 1 ^ 2 ) equation, we find 

Ei\_f 2i(k B d z E + +u> B ( 7B + l)d T E + ) \ AkBI - UBt) 



L °\Pi) \l^d T E + + lB vE_+Z 1 %cl>(\E + \z + 2\E_\z)E, ''' 

2* (-* B Z 1?_ + c^fo + l)d T E_) \ i(kBZ+UBt) 

\^B dT E_ + lB vE + +3 7 |0(|^| 2 + 2\E+\ 2 )E_ 1 



( 

( 



e i(k B z-3Lj B t) 



i(k B z+Zu> B t) 



(3fc B z-w B t) 



(58) 

Of the terms on the right hand side of (58), only the first two are potentially 
resonant, and may therefore give rise to secular growth in time, t. The nonres- 
onance condition required to remove such resonances can be expressed as the 
requirement that the vector coefficients of e'lfe-ust) an( j e -i(k B z+u> B t) \[ e 
in the column space of £o(ui B , ±k B ); see (51). Equivalently, we require that the 
inner product of each of these vectors with the vector (—£02,1, An. i) be equal 
to zero. This yields the Nonlinear Coupled Mode Equations (NLCME): 

t(d T E+ +v g dzE+) + kE_ + r{\E+\ 2 + 2\E_\ 2 )E+ =0 ; (59a) 



|2 1 ol c 1 |2\ 



i{d T E^ ~v g d z E^) + KE + + r{\E-\ 2 + 2\E + \ 2 )E_ =0 . (59b) 



Here 



v g = u'{k B ) = fc'(^)- 1 = "° ^ (60) 

u;b(^-(1 + 7b)) 

is the group velocity, and the coupling and nonlinearity parameters are 

«= , ^ (n2 r 1} ^ r= ^i (61) 
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Our proof of validity of NLCME on time scales of order e 1 requires that we 
solve explicitly for E and P through order e. We solve (58) and obtain 

(*)- E 

V 17 o=±l,±3 V/ l 7 
6=1,3 

such that for (a, b) ^ (±1, 1), 

£ »(|»)) = (sM) ei(0te " M ' (63) 

where S(° ,(> ) is determined in equation (58), and for (a, b) = (±1,1), the right 
hand side is determined by using (59) to eliminate 8tE± from the first two terms 
of (58). Once this is done, these terms take the form of the non-null eigenvectors 
of Co and solving this part of the equation becomes a trivial linear algebra 
problem. In this way we may represent the approximate solution using only E± 
and their Z-derivatives, so that I? and H 8 estimates on solutions to the NLCME 
suffice for proof of the main theorem. 



5. The initial value problem for AMLE and NLCME 

Our proof of the validity of NLCME as an approximation to AMLE requires 
some a priori knowledge of the solutions of these equations. In this section we 
outline the theory of the initial value problems for AMLE and NLCME and 
collect the necessary facts for the proof of the main theorem. 

Both AMLE and NLCME are semilinear hyperbolic systems whose initial 
value problems can be expressed in the form: 

d t <P(t) = -iA$(t) + J[$(t)] , 
<P(t = 0) = -P . ' 

Here, A is a self-adjoint operator on a Hilbert space H and J is a nonlinear 
mapping from H to itself and <Po e Ti. 

We first indicate how AMLE and NLCME can be expressed in this form and 
then show how the general theory and energy estimates can be used to conclude 
the existence of global in time solutions. 

AMLE: 

To write the AMLE system, (6), as a first order system we use the variables: 
E, B, P and Q = d t P. The AMLE system then becomes: 

8 t E = d z B-Q ■ 
dtB = 8 z E; 

8 t P = Q; (65) 
d t Q = - 2evcos(2k B z))P-ujv 1 g(P 7 z)+uj i (n 2 - l)E . 
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We now write (65) in a more compact form. Let 



u = 



B 
P 

W 



and M = 



/ 8 Z -1\ 
d z 
1 
V^ 2 (n 2 - 1) -uit ) 



Then, the full system may be written 

d t u = M.u + uj^e^ {2ev cos(2k B z)P — g(P, z j) , 

where e 4 = (0, 0, 0, 1) T . Thus, 

<P = u, A = iM , 
J[<2>] = w 2 e 4 (2ev cos(2k B z)P - g(P, z)) . 



(66) 

(67) 

(68) 
(69) 



NLCME: 

NLCME can be written in the form (64) with the definitions: 



$ = 



E- 

Of 
10, 

A = —ivgdz — na 1 ; 

2\E-\ 2 )E, 
2|£+| 2 )£- 



jm=<r ((K 



(70) 

(71) 
(72) 
(73) 



We now formulate the general initial value problem (64) as an equivalent 
integral equation: 



$(t) = e- iAt $o + I < 



-iA(t-s) 



J[$(s)]ds 



(74) 



It is elementary to show [38] using the contraction mapping principle that in 
both examples, for any initial condition, $o, in the Hilbcrt space Ti = H 1 , there 
is a maximal time, T max — T max (| \(p Wh 1 ) > 0, and a solution $(t) of (74), which 
is defined for t € [0, T max ), the maximal time interval of existence. The solution 
<P{t) G H 1 for each t G [0,T max ) and the function t ||^(f)|| ff i is continuous 
for t e [0, T max ). Finally, either T max < oo or T max = oo. If T ma;E < oo, then 



Jim \\-P(t)\\m 

t/ 1 max 



(75) 



and we say that the solution <P(t) blows up at time T max in if 1 . As we have seen 
in Theorem 2 of Sect. 3, in the absence of material dispersion, solutions of the 
AMLE system do develop singularities in their gradients in finite time. We claim 
that for both dispersive systems AMLE and NLCME no singularities form: 
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Theorem 3 For initial data in H 1 , T max = oo. That is, AMLE (under Hypoth- 
esis 1 below on the nonlinearity) and NLCME have global in time H 1 solutions. 

To prove this theorem, it suffices to show that if T\ is an arbitrary time, then 
the H 1 norm of any of the components of $(t) satisfies an estimate: 

\\*j(t)\\ H1 <C(7\) . (76) 

The constant, C(T\) may depend on and even grow with T\, but must be finite for 
finite values of T\. To prove (76) we use a combination of the conservation laws 
associated with AMLE and NLCME as well as direct a priori estimates on the 
evolution equations. We consider the cases of AMLE and NLCME individually. 

Proof of Theorem 3 for AMLE: 

We use the formulation for AMLE given in (67) or equivalently (65). Our proof 
makes use of the following technical assumption on the nonlinear term which 
ensures the existence of global solutions for arbitrary size H 1 data: 

Hypothesis 1 There exists a constant C, such that for all z, 

\g(P,z)\ + \d z g(P,z)\ <C\P\ , \d P g(P,z)\<C . (77) 



The first step is to derive an energy estimate for AMLE. Taking the dot 
product of (67) with the vector (E,B, (n 2 - l)' 1 P, Uq 2 (n 2 - 1)~ 1 Q) yields: 



= / z^BzWQdz - -^-j j g(P,z)Qdz (78) 

<C j(P 2 +Q 2 )dz . 
The previous inequality follows from Hypothesis 1. It follows that 

\\m\h<\\to\\h+Ci f \\u(s)\\h ds (79) 
Jo 

for some positive constant Ci and therefore by GronwalPs inequality: 

\\u{t)\\ L 2 < \\u \\ L 2 e Clt . (80) 

Estimates for the L 2 norm of d z il are obtained in a similar manner. We first 
differentiate equation (67) for u with respect to z, and then take the dot product 
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with (d z E,d z B, (n 2 - l^d.P, lu 2 (u 2 - l) _1 a z Q) and obtain: 



d_i r 

dt2J 



E 2 z +B 2 z + ^—P 2 + I Qi | <h 



ev 



(81) 



- y [ 2cos(2k B z)P z Q z - 4:k B sin(2k B z)PQ z } dz 
- -J— I J d P g(P,z)P z Q z dz - J d z g(P,z)Q z dz 

<C'J (P 2 +Q 2 + P 2 + Ql)dz . 

This, together with the above L 2 energy estimate can be used to conclude, by 
Gronwall's inequality 

\\u\\ H i < ||tr || H i e C2t . (82) 

Since the H 1 norm of u grows at worst exponentially, we conclude that T max — 
oo. This completes the proof of H 1 existence for solutions to AMLE. 

Proof of Theorem 3 for NLCME: 

Proposition 1 Let £ = (E + ,E_) satisfy system (59) with initial conditions 
£{Q) e H s for s > l. 9 Then there exists C s = C s (\\£(0)\\ H ,,T) such that 

\\£(T)\\ H s < C s (||£(0)|| hs ,t). Moreover, C{x 1 ,x 2 ) -» as x x -» . 
Proof 

It is easy to see that system (59) preserves the L 2 norm. To obtain this and higher 
LP bounds on E±, we multiply both sides of (59a) by \E+\ 2a E* + and (59b) by 
\E-\ 2a E^, add them, and take the imaginary part, yielding: 

-L^(d T (\E_\ 2 °+ 2 + \E + \ 2 °+ 2 )) +vd z (\E + \ 2 °+ 2 - \E_\ 2 °+ 2 ) 

+ ir{E+E*_ + E-E* + )(\E+\ 2a - |£L| 2ct ) = . (83) 

If a = 0, then the last term is identically zero, showing that \\£ || 2 is conserved. 10 
If a > 0, then we may bound \\£\\ using Gronwall's inequality 

| 2CT+2 < ll£lL _,=(*+i)T ■ (84) 



\\S\QZ < Hfo|| 2ff+2e c(CT+1)T ; 

H^ll2cr+2 — l|£o|| 2(T +2 C 



9 In our proof of Theorem 1, we use this result for s < 3. 
10 Recall that \\£\\ P p = f (\E+\p + \E_\P)dZ. 
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Letting p = 2a + 2, this is just 

\\£% < \\£v\\ p e cT/2 • (85) 

As c is independent of p, this estimate holds for L°° . 

The L°° bound can then be used to bound growth rates of the L p norms 
of d z E± in terms of T. Taking ^-derivatives of the NLCME, and performing a 
similar calculation with a — yields: 

^\\d z et<c\\e-d z it 

<c\\£\\ 2 LOO \\d z £\\ 2 L2 ( 86 ) 

<c\\£^ L ^ T \\d z £\^ , 

so that 

ll^ltfi < ||3>||*ie c < eCT - 1 > • (87) 

This shows that we can indeed bound the solutions of NLCME in H 1 , and 
control them for times T = 0(1), i.e. t = O(-). Proceeding similarly, we can 

ce cT 

derive bounds in higher Sobolev spaces, specifically H 2 norms like e ce , and 

H 6 bounds like e ce , thus completing the proof of Prop. 1, and hence, by 
the comments preceding Theorem 3, of that theorem. 

6. Validity of NLCME for times, t, of e>(e -1 ); proof of Theorem 1 

We shall work with the formulation of AMLE given in (65). 

We proceed under the following hypothesis concerning the nonlinearity g(P, z) 
and its derivative dpg(P, z) for small P: 

Hypothesis 2 Assume g has partial derivative of order < 4 with respect to P 
and has one partial derivative with respect to z which is continuous. Assume 
further thatg{0,z) = (d P g)(0,z) = {d P g){0,z) = 0, and <f> = -^(d%g)(0,z) + 
and is independent of z, and make the analogous assumptions for d z g Therefore, 
there exists a positive constant, C, such that for all z and all P\,Pi with \P\\ + 
I-P2 1 sufficiently small: 

| 5 (Px + P 2 ,z) g{P u z)\ < C (|P!| 2 + |P 2 | 2 ) |P 2 | ; 
\d P g(P 1 +P 2 ,z)-dpg(P 1 ,z)\<C(\P 1 \ + \P 2 \)\P 2 \ ; (88) 
\d z g(Pi + P 2 ,z) - d z g(P u z)\ < C (|Pi| 2 + |P 2 | 2 ) |P 2 | . 
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To obtain an approximate solution of (65), we require, in addition to our 
approximations of E and P, approximations to B and Q through first order 
in e. We use the relation d t B(t, z) = d z E(t, z) to obtain the relations 



dtBo = d z Eo ; 

d t B x = -d T B + d z E 1 + d z E . 
Also, using Q — d t P, we find 

Qo = d t Po ; 

Oi = dtPi + d T Po ■ 

We may then define 

Xl m = e2 (X + eXj) for X = E,B, P, or Q 
and write our approximate solution to (65) as: 



/E £ \ 

app \ 
rye 

app 
pe 

app 

\QappJ 



The full solution to AMLE may therefore be written as 

u = if app {t,z;T,Z) + eRc(t,z) , 

where 



(89) 
(90) 



(91) 
(92) 



(93) 



(94) 



(R%\ 

R £ P 



(95) 



(96) 



denotes the error term. To prove the main theorem it suffices to prove that for 
any T > 0, R e remains bounded of order one, in an appropriate norm, uniformly 
for e sufficiently small and < t < T /e. 

We now derive the equation for R e . Viewing i, z,T, Z as independent vari- 
ables, (67), the equation for u £ can be rewritten as: 



d t u — Mu - eNu + cjQe 4 (2ev cos{2kBz)P — g(P, z)) , 



where 



I d T -d z \ 
-d z d T 
d T 
\ d T ) 



(97) 



(98) 
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To obtain an evolution equation for R e , we substitute (95) into (97) to obtain 

d t {ul pp + e&) = M(u £ app + eP>) + eN{ul pp + e&) 

+ wge 4 (2e^cos(2fc B z)(P a £ pp + eR P ) - g{P e app + eR P , z)) . 

(99) 

We may then eliminate from this two equations obtained during the multiple 
scales expansion: 

(d t - M)u = ; (100a) 
(d t - M)u! = - Afu + 2u 2 e 4 cos(2k B z)P + Lu 2 e 4 (f>P ( 3 (100b) 

to leave an equation for the evolution of the error alone 



(d t - M)R £ = - eiNui + 2ve 4 cos(2k B z)(ei Pi + eRp) 
+ Lufai-e-'giP^ + eRp, z) + £ ^P 3 ) 



(101) 



To this we add and subtract e ujQe 4 g(P £ z) to obtain 



app 

1, ,2; 



(d t - M)R e = -eiAfui + 2ve 4 cos(2fc B z)(e^Pi + eRp) 
+ £- 1 ^e 4 ( ff (P a e pp , z ) - g{P e app + eRp, z)) 
+ c 2 e 4 (- £ - 1 .g(P a £ pp , z) + £ ^P 3 ) (102) 
= 2e^e4 cos(2fcsz)Pp 

+ e" 1 ^e 4 (.9(Pap P , z) - g{Pa PP + eRp, z)) + e~ x p , 

where 

p = -ei/Vtfi + e*wg2i/cos (2fc B z)P 1 + (efyf§ - .g(P a £ pp , «)) (103) 

is the residual, essentially the amount by which u e app fails to solve (67). 

We now consider the formal size of the second and third terms on the right 
hand side of (102). Since P a £ pp + eRp = e? (p + ePi + e? Rp) and g(P, z) ~ 



</P 3 , the second term of (102) is 0(e). We further note that since u e app is an 
approximate solution through order e, the residual p is formally of order e 2 and 
so the third term of (102) is 0(e). It is in order to control this final term that 
requires us to calculate the approximate solution including terms formally of 
order £2 and also to require that E± be in the Sobolev space H 3 . 

From this discussion, we expect that for times of order e^ 1 , R £ will be bounded. 
For convenience, we introduce a notation which makes explicit the expected size 
of the residual: 

ef=s~ 1 p. (104) 
It is then clear that in order to bound R £ we will first need to bound r. 
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Proposition 2 (Estimation of the residual,) Let (E + ,E_) be a solution of 
the NLCME system (59). Then there exists a constant c > depending on k and 
the Sobolev norms of E± of order up to three, but independent of e such that for 
allO<t< T e-\ 

M\ L oc <ce x l 2 ; (105) 
\\r\\ Hl < c . (106) 



This proposition is a simple consequence of the explicit expression for f (defined 
in terms of p) given in (103), and of Prop. 1. 

Proposition 3 Let (E + ,E_) be a solution of the NLCME system (59). Then 
there exist e > and Co > s.t. if < s < Eq, the solution of (102) satisfies 
\\R £ \\ H i < Co for all0<t< Toe" 1 . 

These propositions imply Theorem 1. 



6.1. Proof of Proposition 3: Estimates on the error R £ 
Recall that the evolution equation for R e is given by: 

d t R £ = MR £ + e 4 uj 2 ( y 2evco&{2k B z)R e P - J (g{P, z) - g(P^ pp , zjfj + er , 

(107) 

Motivated by the energy estimates used in Sect. 5 introduce the weighted 
norms 

//■oo r>s 2 p^s 2 

if* • A& dz = / (i?| 2 + R% 2 + + 9 ) dz ■ 

J-oo n 2 -l uj£(n 2 -l) J 

(108) 

\\\R £ \\\m = lll^lll 2 + lll«lll 2 , (109) 

where the weight A is given by the matrix: 

1 1 



A = diag 1, 1, 



Vu 2 (n 2 - 1) 



These norms are clearly equivalent to the standard L 2 and H 1 norms. 
Then we have 

l^ ll|i?£|!|2 J ( 2 ^ C0 ^ k B-AR% + \[9{P'a PP ^)~9{Pa PP + eR e P ,z)] )R 

+ e j R £ ■ Ar dz. 

(110) 
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Similarly, differentiation of (107) with respect to z and left multiplication by 
(d z W) T A yields: 

^Hl^lll 2 = ^ZJ J {2cos(2k B z)d z R e P d z R £ Q - 4k B sm(2k B z)Rpd z R s Q ) dz 
(dpgiP^pp, z) - d P g(P^ p + eR P , z))d z P s app d z R s Q dz 
dpg(P ap p + eRp, z)d z RpR e Q dz 
+ -J (d z g(P! pp , z) - d z g(P s app + sR P , z))d z R £ Q dz 



1 

£ 



7 



+ e / d z R e ■ Ad z fdz . 



(Ill) 

Application of Hypothesis 2, the L°° bound on solutions of NLCME of Propo- 
sition 2, Sobolev's inequality 11 [22] and interpolation yields 

^lll^ £ IH 2 < Cie|||fi e ||| 2 + C 2 e 2 |||i? £ |||^i +e\\r\\ L2 \\\&\\\ . (112) 

and 

^Ill^^-IM 2 < | 1 1 ^ C^^ 2 1 1 1 1 1 J^^ <5»^^|| \\\R e \\\m ■ (113) 

Estimates (112) and (113) imply 

j t \\\R e \\\m <C (e\\\R £ \\\ 2 m + e^R^ + e\\r\\ Hl |||i£|M • (114) 

If R s (0) = 0, then by equation (107), ^ for t > 0. We therefore 

assume ||i? e (0)|| > 0. We may then let = |||i? £ |||(i). Then (76) and the H 1 
bound on r from Prop. 2 implies 

§<C £ (l + C + eC 3 ) . (115) 
This differential inequality is easily solved and we conclude that for any T > 0, 
\\R E (t)\\m < C |||i? e (i)lll < C(T ) for < t < T^ 1 . (116) 
Finally, we note that, as 

u = unlcme + +eR £ , (117) 

then 

3 — * 

\\U — UnlCME\\ fji < £ 2 \\ui\\ H i +e\\R 6 \\ Hl . (118) 

The estimates on E± guarantee that the first term on the right hand side is 0(e) 
and Prop. 3 guarantees that the second term is 0(e). Thus 

\\u-u NLCM e\\ h i < Ce for < t < TqE^ 1 . (119) 

This completes the proof. 

11 ll/llloc <qi/!l L 2||a z /|| L2 
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7. Numerical Demonstration gap soliton propagation and decay 

We initialize a wavepacket for AMLE with many oscillations and an envelope 
whose form is constructed using the gap soliton solution to the NLCME. The gap 
soliton decays exponentially away from its 'center'. We perform the simulations 
with periodic geometry. The period is chosen to be several gap soliton widths 
so that the solution is well localized away from the artificial period ends. The 
gap soliton initial condition is initialized within the central region of the domain 
so the localized structure is essentially unaffected by the boundary and propa- 
gates as though it were on an infinite spatial domain. In addition, we compute 
the evolution of a solution to the NLCME with corresponding envelope initial 
conditions, and use the formulae (54), (89), and (91) to construct approximate 
solutions to AMLE for comparison. 

7.1. Numerical Methods 

We use a "method of lines" approach, meaning that we first discretize in the spa- 
tial dimension, yielding a set of ordinary differential equations in t for the values 
of the solution at the discretization points. We compute solutions to AMLE 
as the first order system given in (65). We restrict our computational domain 
to a finite period interval and discretize with about 16 points per wavelength. 
Derivatives arc computed spectrally using discrete Fourier transforms [20] . That 
is, suppose T is the discrete Fourier transform, £ is the dual variable to z, and / 
is a vector of discrete values of /. Then the approximate derivative is given by 

d z f = (T-^T) f . 

The spatial discretization of system (65) may now be treated numerically as 
a system of ODE's. A fixed-step fourth order explicit Runge-Kutta method is 
used to integrate the resulting system in time. Recall that an n-stage explicit 
Runge-Kutta method for a general evolution equation 

y(t) = f(y,t) 

is given by [21]: 

fc i = f{yk,tk) ; 

i-1 

fcj = f(yk + aiAt, tk + At bijkj) for i = 2,n ; 

n 

Vk+i =Vk + ^2 cik t . 
i=i 

Explicit methods tend to impose stability restrictions on the allowable step 
size for the time integration. However, in the case of the AMLE (65), this is 
simply that At < CAz, which is a very mild restriction compared, for example, 
to the heat equation, for which At < CAz 2 . For most of our simulations we work 
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with about 20 points per wavelength, which gives Az 



\ and a comparable value 



for At. 



Empirical convergence tests show the method to have fourth-order conver- 
gence in time, and constants of motion are also computed numerically and are 
shown to be conserved to 8 or 10 digits. A similar method is used for NLCME, 
though the accuracy is far less crucial as the solutions contain far fewer oscilla- 
tions and vary on a slower time scale. 

7.2. Numerical Verification 

To numerically verify and explore the limits of Theorem 1, we solve both AMLE 
and NLCME under the SVEA scaling and compare u with y/euo by monitoring 
the quantity: 



We do this for two values of e, and check that the agreement scales appropriately 
as e is decreased. This is done for £\ = ^ and Si = so the error should be 
reduced by half between the two runs. For the purposes of verification, we take 
much larger values of the refractive index contrast e than would be used in a 
physical experiment. 

In Figure 7, we show a typical initial condition for the electric field E. The 
parameters used for this and all the numerical experiments described in this sec- 
tion are given in Appendix C. The envelope in this figure is 256 wavelengths long, 
and it is generated from a simulation with e — 1/64. The shape of the electric 
field envelope as the wave propagates clearly illustrates the effect of the periodic 
medium on propagation. In Figure 8 we sec that the electric field envelope (com- 
puted from "full" AMLE solutions) is "two humped" with the amplitude moving 
forward and backward between the humps at a faster rate than the envelope it- 
self moves forward. In the same figure, we plot the location of the maximum of 
the electric field, and it becomes clear that the electric field maximum moves 
forward unsteadily, interrupted by a sequence of backward jumps. Also plotted 
in this figure is the location of the energy density maximum, 12 which propagates 
more smoothly, since the contribution from the different fields is averaged. 

Figure 9 shows the location of the electric field envelope at the beginning, 
middle, and end of the computed evolution period. This figure shows both the 
envelope computed from the AMLE and also the approximate envelope com- 
puted using the NLCME. To the "eyeball metric", the agreement appears to be 
quite close. More quantitatively, the success of this procedure is measured by an 
error-scaling factor given, for any norm, 



Then the numerics verify the asymptotic procedure if the scaling factor is equal 
to one. Figures 10 and 11 show that computed in L 2 the error scales in agreement 

12 The energy density is the integrand on the left hand side of equation (78). 



Error e (i) = u £ — \fs~Uo . 



(120) 



Error Scaling Factor = log 2 



\\Error £ \\{T e ) 



\\Errors\\(Te) ' 
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Fig. 7. The initial data for the electric field satisfying the SVEA. 
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Fig. 8. (Left) The motion of the electric field by reflection and nonlinear regrouping. Between 
(a) and (c) the envelope moves forward, between (c) and (e) it is reflected backwards, and at 
(f) it has begun propagating forward again. (Right) The location of the maximum of E (solid 
line) and the location of the energy density maximum (dashed line) as a function of time, 
showing the effect of reflection off the grating (computed using full AMLE system). 



with Theorem 1, but that the L°° error is reduced by a factor of 22 . A general 
scaling argument shows this is reasonable. Consider a function f(z) and let 
f £ (z) = ei.f{ez), then ||/ £ || 2 = e||/|| 2 , while WfeW^ = e'H/IL- It appears that 
using H 1 estimates to control L°° estimates has cost us half a power of e in our 
approximation of the error R e . 



7.3. Very Long Time Behavior 

The error estimates of Theorem 1 tell us that the solution constructed from the 
NLCME and the full solution to AMLE should agree for times on the order of 
£ _1 . As a practical estimate, this may cause us some concern, as the width of 
the solitary wave is also 0(e~ 1 ) so that on these time scales, the distance of 
propagation is the same order of magnitude as the width of the solution. It is 
therefore of interest to run our simulations for long times, to see if the NLCME 
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Electric Field Envelope 
t=0 




-1500 -1000 -500 500 1000 1500 



Fig. 9. The envelope of the electric field at the beginning, middle, and end of the computed 
evolution. Computations of both the AMLE envelope and its NLCME approximation are 
shown. 



continues to provide a good approximation beyond what we have proven, or if 
the approximation breaks down completely. 

We run the simulation with e = 1/32, and with w = 4, allowing the evo- 
lution to continue to t = 12000, which is certainly larger than 0(e~ 1 ). By this 
time, the L 2 norm of the error is of similar to the L 2 norm of the field itself, and 
has stopped growing. Figure 12(a) shows that the envelopes of the full solution 
and the approximation no longer agree, but that they lie in approximately the 
same location. In Figure 12(b), we see a blow up of the electric field and its 
approximation via the NLCME which shows that the two solutions arc com- 
pletely out of phase with each other, so that pointwise estimates will not show 
any agreement between the solution and the approximation. Figure 12(c), how- 
ever, shows that the energy density of the full solution and the approximation 
continue to match very well. The solution has propagated about twenty times 
its own width (full width at half maximum or FWHM) , and the centers of the 
two energy density plots are separated by about one fifth of a FWHM. This is 
encouraging, as it suggests that, although the estimates of Theorem 1 no longer 
hold, the approximation and the full solution have basically remained together. 

Although the description via the NLCME has broken down in the above dis- 
cussion, the electric field has maintained the basic structure of a slowly modu- 
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Fig. 10. Scaling of the error in L 2 as a function of the scaled time. 



lated plane wave. Eventually, this very structure will break down and the solitary 
wave may itself break apart. This is shown in Figure 13 where a solitary wave, 
moving to the right, steepens at its trailing end and then begins to break up, 
while falling behind the AMLE envelope. In this figure, the parameters are as 
in Figure 12, except that ujq = 1, this has the effect of decreasing the number of 
oscillations contained within the FWHM of the envelope from 60 to 10, so that 
the separation of scales is much less pronounced. This much narrower envelope 
breaks up much faster than the solutions shown in previous plots. Although we 
have no precise measurement of this breakup time, it appears to happen on a 
time scale t ~ 0{s~ 2 ). 

8. Summary and Discussion 

In this paper we considered the propagation of high intensity light through a one- 
dimensional periodic structure. This was modeled by the anharmonic Maxwell 
Lorentz equations (AMLE), which incorporate the effects of material dispersion 
due to finite time response of the polarization field, photonic band dispersion 
due to the periodic structure, and nonlinearity (intensity effects). We first gave 
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L°° error 




Fig. 11. Scaling of the error in L°° as a function of the scaled time. 



a detailed discussion of how these effects act individually and in various combi- 
nations, also providing some numerical illustrations. We next considered AMLE 
solutions with spatially localized initial conditions, nearly monochromatic at the 
Bragg resonant carrier frequency and such that the effects of dispersion and non- 
linearity are balanced. We proved that, over time scales of interest, the backward 
and forward propagating field envelopes satisfy nonlinear coupled mode equa- 
tions (NLCME), which we derived from AMLE using multiple scale analysis. 
We also derived rigorous bounds on the deviation of the NLCME solutions from 
those of the original Maxwell-Lorentz model. Theorem 1, which describes this, 
is the main mathematical result of the paper. We demonstrated its validity and 
probed its limitations via numerical simulation, as well as verifying that the 
ordering assumptions assumed in our analysis are consistent with the physical 
parameter magnitudes characteristic of experimental studies. 

Two directions of great interest are the the study of nonlinear phenomena 
in multidimensional photonic structures [3] and the extension of the present 
analysis to the case of more general inhomogeneous structures. The multiple scale 
techniques and the analysis used to obtain Theorem 1 can be applied to more 
general structures, e.g. periodic structures defined by a general Fourier series, 
index variations which are slow modulations of those considered, and "deep 
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Fig. 12. The full (solid) and approximate (dashed) solution at t = 12000 for (a) The electric 
field envelope, (b) the full electric field, blown up from the box in (a), showing that the two 
solutions are out of phase, and (c) the energy density which agrees quite well. 



gratings" . In the case of deep gratings, where the variation of the refractive index 
is not small, this requires the use of a multiscale expansion ansatz describing the 
slow modulation of Floquct-Bloch waves [10], rather than the plane waves we 
have used in the case of a system which is nearly translation invariant in z. 

A number of issues arose in our study we which presently discuss and raise 
as questions meriting further investigation: 

(1) Numerical simulations suggest that NLCME continues to acceptably predict 
the location of the field energy on time scales for which the estimates used in 
proving Theorem 1 break down. Fig. 12 shows that the coupled mode theory fails 
to predict the location of the individual peaks of the carrier wave while continuing 
to predict the location of the energy. It would be of interest to investigate whether 
there exists a weaker, more general framework, in which the AMLE solution is 
well described by the NLCME solution. 

(2) The very long time simulations described in Sect. 7 indicate a degradation 
of the gap soliton due to wave steepening and the radiation of energy away from 
the soliton core. It would be of interest to derive higher order model equations 
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Fig. 13. The envelope of the solution with very few oscillations steepening and then breaking 
up 



which describe these phenomena and agree with the full solutions to AMLE on 
longer time scales. 

While it is possible to find longer-time envelope equations by starting with 
smaller (O(e)) initial conditions and introducing a third time scale T 2 = e 2 t 1 our 
primary interest is to investigate the validity of the NLCME system already in 
wide use by experimentalists. Rigorous results for such longer-time systems in 
other contexts are given in [23, 31]. 

(3) For the one-dimensional non-dispersive model with nonlinearity, we have 
seen that wave steepening and shock formation occurs. This situation appears 
to persist in the presence of periodic structure. It would be of interest to extend 
the Lax-Klainerman-Majda theory [28] described in Sect. 3.4 to include the case 
of equations with periodic or more general inhomogeneous variable coefficient 
terms. As noted, for (localized) SVEA initial conditions with carrier frequency 
in Bragg resonance with the medium, photonic band dispersion is significantly 
stronger than material dispersion and the gap solitons arise due to a balance 
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between the former and the Kerr nonlinearity (indeed, the NLCME may be ob- 
tained as a Galerkin truncation of the infinite system of equations (39)) derived 
in this case. Inclusion of material dispersion (wo < oo in AMLE) regularizes 
shocks; see Theorem 3. Are there subtle regularizing effects provided by pho- 
tonic band dispersion alone? 

(4) In the full three-dimensional waveguide problem, one must also take into 
account waveguide/mode dispersion and polarization mode dispersion. In this 
case there is an interplay between the mechanisms of diffractive spreading (regu- 
larizing), geometric confinement of the field (tending to onc-dimcnsionalize and 
therefore singularize the propagation) , modal dispersion (which takes higher har- 
monics off resonance and therefore possibly regularizes) and nonlinearity. The 
interplay of all these effects remains unclear. It would be interesting to extend 
the results of [12, 23, 31, 40, 39] to situations with periodicity and nontrivial 
transverse geometry. 

A. Dimensionless Quantities 

In this appendix we nondimcnsionalizc AMLE and isolate the key nondimen- 
sional parameters. We then define dispersion lengths and nonlinear length whose 
balance specifies the conditions under which a soliton is expected to form. Fi- 
nally, using the experimental parameters of Egglcton ct. al. [15], we calculate 
our dimensionless quantities and verify the applicability of AMLE. 

We begin with the AMLE system written using dimensional variables and de- 
rive a nondimcnsional version of AMLE. We then use physical parameter values 
gleaned from the literature in order to find approximate sizes of the nondimen- 
sional parameters. Primed variables represent nondimcnsional quantities and 
unprimed variables dimensional ones. We use the standard notation [X] to rep- 
resent the units of X so that X = [X]X' for any variable X. 

The AMLE written in dimensional variables are: 



H d t D = d z B; d t B = d z E ; 
D = e a E + P ; 

tio 2 d? p + (! - 2An cos(2k B z))P - 4>P 3 = e aX 



(121a) 
(121b) 
(121c) 



We begin, as usual, by eliminating the magnetic field B to obtain 



(122) 



We now introduce nondimensional (primed) variables: 



t = Tt'; z = Zz' ; 
E = £E' ; P = e a £P' ; D = e £D' ; 

r. k B - ujq 



(123a) 
(123b) 
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where the calligraphic letters represent dimensional magnitudes. To explicitly 
display the expected scaling, we write 

An = ev and 4> = £< t , (124) 
(eoc-J 

where v and (f>, along with x^\ are dimensionless and 0(1), and the fields in 
(123b) are also all 0(1). 

Substituting these new variables into equations (122), (121b), and (121c) and 
eliminating common factors yields: 

^,D' = ±dlE> ; (125a) 
D' = E' + P' ; (125b) 
-\d 2 ,P' + [1 + 2ev cos(2fc B z')]P' - ecfiP' 3 = X (1) E' . (125c) 



LUq 2 



Letting 



T=|, (126) 



we have that (125a) becomes 

d 2 tl D' = d 2 zl E' . (127) 

The system (127), (125b), and (125c) comprise the dimensionless AMLE system; 
see also (6). 

A.l. The material frequency Cjq and the electric susceptibility 

At low intensities the relation between P and E is given by the Lorentz model: 

-Lp tt + P = e QX {1) E . (128) 

In the time-frequency domain this implies: 

P(w) = eo ^" 2 E(u) , (129) 

where f{uo) = J e-' iut f{t) dt. 

In a general linear setting we have, 

P(u)= X (1) (u ) )E , (130) 

where the (frequency dependent) index of refraction, n(uj), is related to by 
the relation: 

n 2 (uj) = l + x (1) H • (131) 
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A standard model for x^ ( w ) m the optics literature [2] is the Sellmeier model, 
which approximates x^( w ) by a function of the form: 

N 2 (1) 

X (1, H-oE^ , (132) 

»=i 1 

where u>i, the model resonant frequencies of the medium and x'P are determined 
by a data fit. For silica glass, a good fit with experimental data is found with 
N = 3. The Lorentz model corresponds to N = 1, so we take the term in the 
N = 3 expansion corresponding to that frequency, u>i which is closest to the 
input carrier frequency. Below, we use this to determine the values of u>o and 
X^ in the Lorentz model. 



A.2. The Electric Field Strength £ 

Most optical physics literature reports field strength in terms of the intensity, /. 
The electric field strength is given in terms of the intensity by [2] : 

£ 2 = — . (133) 
e cn 

where n is the (nondimensional) refractive index, related to the linear suscepti- 
bility, by (131). 



A. 3. The coefficient of nonlinearity, <fi. 

We consider the instantaneous limit of the basic equation, with no grating, i.e. 
An = 0: 

P-4>P 3 = e x (1) E . (134) 
We may invert the above relation for small E and write 

P = e ( X (1) £ + X (3) £ 3 + ...) , (135) 

where 

,0) 



e 2 X^ 3 



(136) 



Then the nondimensional quantity is given by 



£<t, = W)f- (137) 

The third order susceptibility x^ is related to the nonlinear refractive index, 
n 2 or n ! 2 by the relation ([2], page 40, equation (2.3.13) and page 582, equation 
(B.2)): 

^ (3) = 8mi2 = Ae cn 2 nl, 
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Finally, since I = \t§cnE 2 ([2], page 582, equation (B.l)), we have 

Mnn{ (139) 



Y 3( X «)3 ' 
A.Jf.. Parameter values of physical experiments 

To form the anharmonic oscillator equation for the polarization, we need four 
constants: the susceptibility the nondimensional frequency luo, the index 
modulation An, and the cubic coefficient <f>. 

The Susceptibility \^ an( l the nondimensional frequency lu 

First we must find the characteristic time scale T. Typical experiments are 
performed using laser light with wavelength of approximately one micron. We 
define the characteristic length and time so that ks ~ 1, but for convenience in 
the paper refer to ks- Accordingly, we take 

i x in -6 

Z = m w 1.6 x l(T 7 m ; (140) 

T = — w 5.3 x l(T 16 s . (141) 

c 

Next, we must find the dimensional frequency of the oscillator. For silicon 
glass, one has ([2], page 7) 

Co {) = 1.6 x lC^V 1 ; (142) 

X {1) = .41 . (143) 

The nondimensional resonant frequency is then given by 

luq = u T w 8.6 . (144) 

The index modulation An = ev 

Egglcton et. al. [15] give an approximate value of 

An w 3 x 1CT 4 . (145) 

The nondimensional nonlinearity coefficient, (j> For this we need the in- 
tensity, which in [13] is given by 

Iw2xl0 14 W/m 2 (146) 

and the nonlinear refractive index ([2], pages 582-583): 

n{ = 2.5 x 1CT 20 m 2 /W . (147) 

The linear refractive index is obtained from (131) and (143): 



n w 1.2 



(148) 
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From (139) we have 

e<p w 2 x 10~ 4 . (149) 
Therefore, by choosing the small parameter 

e = 1(T 4 , (150) 

we arrive at 

0w2 (151) 

and 

i/w3 . (152) 

Therefore the approximate nondimensional polarization equation (125c) may be 
written: 

C'(10- 2 )9 t 2 ,P'+[l + C'(10- 4 )cos(2fc B z , )]P' + C , (10- 4 )P' 3 = C'(l)£; , . (153) 

We see that the nonlinearity and the dimensionless grating effectively balance 
each other. This justifies our e-dependent scaling of the dimensionless AMLE 
system and the solution. We note that this scaling assumes that E' , P', and D' 
are 0(1) quantities, see (123b). In the main text, we take E, P, and D to be 
0(y/e), thereby effectively introducing the factor of e multiplying <f) m (125c) 
which is absent from (6c). 

Note also, that while w ~ 2 , the coefficient of df,P', is small, it is roughly 100 
times the grating strength, i.e. eui 2 ~ 10~ 2 . The significance of this can be seen 
as follows. Were we to expand the electric field as in Sect. 4, to all orders in e 
we would have: 

oo 
i=0 

Inspection of the hierarchy of equations for Ei reveals that 

Ei ~ ujo 2t . 

This suggests that E £ is well-approximated by Eq provided euio 2 *C 1. The 
experimental regime discussed satisfies this criterion. 

B. Calculation of the dispersion and nonlinear lengths 

In the design of an experiment to observe gap solitons, the matter of the forma- 
tion length is important. Laser light injected into an optical fiber will have an 
approximately Gaussian profile. One is therefore interested in the distance over 
which one can expect a soliton to form. Solitons are understood to form due to 
a balance of dispersive and nonlinear effects. Dispersion acts by broadening a 
pulse and radiating high frequency components away, while a Kerr (focusing) 
nonlinearity acts to concentrate energy. We presently give a heuristic discussion 
of this balance. 
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Material dispersion length, Zd, material 

Recall that the (material) dispersion relation associated with the finite time 
response of the medium to the field is: 

l-(-) 

The dispersion of a wave packet, with frequency content concentrated in an 
interval of width e about lob is governed by Fourier integrals of the form: 

I(z,t) = fe^W-^) f (^R )duj ; (155) 

where / is a localized function of frequency. Expansion of k(u) about to = uob 
yields 

I{z',t') <~ e ^ B z'~u B t') J e i(uj~LU )(k' (u> B )z' -uj B t') e i k " { Z B) (u-LUnft' - ^B ^ ^ 

= 0{{k"{oj B )e 2 z 1 )- 1 -) , z' = 0(s- 2 ) . 

(156) 

Thus, a localized pulse disperses due to the finite time response of the medium 
over a dimensionless distance z' of order e~ 2 k" (ujb) -1 ■ Noting that k"(uJB) = 
0(uJq 2 ), we have: 

ZD, material = 0(e~ 2 k" (oj b ) _1 ) = 0(w 2 e~ 2 ) wavelengths. 

Using the physical parameter values discussed in Appendix A, we find that 

Z D, material ~ 7 km . 

Photonic band dispersion length, ZD,band 

Linear dispersion due to the periodic structure (photonic band dispersion) is 
governed by the linear coupled mode equations (27) . The dispersion of the wave 
envelope is then expressed in terms of generalized Fourier superpositions of 
Floquet-Bloch waves: 

I(z',t')~ J E(z',t',;k B + eQ)f(Q)dQ 

^ e i(k B z'-0J B t') e -W(0)et' J e iQ(ez'-fi'(0)et') & - (0)Q 2 et' j (Q^g 

= O{{Q"{0)et)-h) , t'^Oie- 1 ) , 

(157) 
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where Q{Q) is given by the dispersion relation (29) which disperses to zero over 
a distance 

2 D ,ba„d = OdeV'iO)}- 1 ) = O(^) wavelengths. 
Physically, this gives 

2D,band ~ 1 CHI 

which is six orders of magnitude shorter than zd, material- 



Nonlinear length, znl 

A measure of the distance, znl, over which nonlinear effects play a role can be 
obtained by considering the coupled mode equations in the absence of dispersion. 
If Eq denotes the electric field amplitude then we have: 

{d T ±v g d z )E = -ir\E\ 2 E (158) 

with solution E = e -^£ 2 (z-v g T) £ = e -ir£ 2 (ez' -v g st') £ _ for somc constant £ 
Therefore, 

^nl = (e^) 1 ~ 1 wavelengths 

which gives 

znl w 1 cm 

which balances the band dispersion length ZD,band- 
Balance of nonlinearity and dispersion 

Note that zd, material is longer than ZD,band by a factor of order e ; for frequencies 
near the band edge, the dispersion due to the periodic structure is much stronger 
than material dispersion. 

Therefore, in order to achieve a balance between dispersive and nonlinear 
effects over a short distance we must equate ZD,band and znl- This gives 

K 1 

— <~ — or K(f) ~ e . (159) 

£ (f> 

By (139) this is gives the intensity ensuring a balance of appearance of nonlinear 
effects within a (photonic band) dispersion length, ZD,band: 

I r^i 

8 nn^K 

which works out to 

/ = O(10 14 ) W/m 2 , 
in line with the experiments described in [13]. 
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C. Computational Details 

C.I. Computations in Sects. 3.4 and 3.5 

All figures in these sections are initialized as a single normal mode of wavenum- 
ber k = 1 of the linearized form the AMLE, 6, with v = and magnitude sfs. We 
should note that the numerical simulations use very large values of e compared 
to the physically appropriate value e — 0(1CP 4 ) derived in Appendix A because 
performing the simulations for AMLE with such small e would be computation- 
ally infcasible. The other parameters are given by 

n 2 - 1 = 1 and <j> = 1 . 

For Fig. 4, we use a frequency of w = 1000 to illustrate the behavior near the 
limit of instantaneous polarization before the onset of a shock. In Figs. 5 and 6, 
we use ojo = 50 and uo = 100, respectively to show the role of dispersion in 
regularizing the shock. 

C.2. Computations in Sect. 7 

All calculations in this sections are performed with the following parameter 
values: 

1 

£ ~ 32 ; 
k B = 1 ; 
n = 1.19 ; 
4>=\ ■ 
v = 1 . 

In Figs. 7, 9, and 12, a material frequency of 

luq = 4 

is used, while in Fig. 13, we use the value 

OJQ = 1 • 

In all calculations, the coefficients of the NLCME are derived from the above 
parameters, and, as initial conditions, we construct a solution from the gap 
soliton with parameters 

v = .9 and S = .9 . 

To create the graphs in Figs. 10 and 11, we also compute the evolutions with 
all parameters as above except with e = wj. 
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